

namespace hms::swe::hllc
{


template<typename RawFlux_out, typename D, typename V>
void setRawFlux(
	RawFlux_out& rawFlux_out,
	const D& d,
	const V& vNormal
){
	rawFlux_out.row(0) = massFlux(d, vNormal);
	rawFlux_out.row(1) = momentumFlux_normal(d, vNormal);
}

template<typename RawQ_out, typename RawFlux, typename D, typename V>
void setRawQ(
	RawQ_out& rawQ_out,
	const RawFlux& rawFlux,
	const D& d
){
	rawQ_out.row(0) = d;
	rawQ_out.row(1) = rawFlux.row(0);
}


template<typename Flux_out, typename F, typename V, typename N>
void setFlux(
	Flux_out& flux_out,
	const F& f,
	const V& vPara,
	const N& normal
){
	static_assert( Flux_out::RowsAtCompileTime == 3 );
	static_assert( F::RowsAtCompileTime == 2 );
	static_assert( N::RowsAtCompileTime == 2 );

	assert( rawFlux.cols() == flux_out.cols() );
	assert( normal.cols()  == flux_out.cols() );
	assert1xN_or_arithmetic( vPara, flux_out.cols() );

	auto n = [&](Index row) -> decltype(auto) {
		if constexpr ( isDenseType<V>() ){
			return normal.row(row).array();
		} else {
			static_assert( N::ColsAtCompileTime == 1 );
			return normal[row];
		}
	};

	flux_out.row(0) = f.row(0);
	flux_out.row(1) = n(0) * f.row(1) - n(1) * f.row(0) * vPara;
	flux_out.row(2) = n(1) * f.row(1) + n(0) * f.row(0) * vPara;

	// auto sign = [&](){
	// 	if constexpr ( isDenseType<F1>() ){
	// 		return Array2s{-1,1}.replicate( 1, normal.cols() );
	// 	} else {
	// 		return Array2s{-1,1};
	// 	}
	// }

	// flux_out.template bottomRows<2>() =
	// 	normal.array() * repl<2>(f2) +
	// 	sign() * normal.array().colwise().reverse() * repl<2>(f1 * vPara);
}

template<EdgeRange orien, typename Flux_out, typename F, typename V>
void setFlux_ortho(
	Flux_out& flux_out,
	const F& f,
	const V& vPara
){
	static_assert( Flux_out::RowsAtCompileTime == 3 );
	static_assert( F::RowsAtCompileTime == 2 );


	using E = EdgeRange;
	if constexpr ( orien == E::horizontal ){
		flux_out.row(0) = f.row(0);
		flux_out.row(1) = - f.row(0) * vPara;
		flux_out.row(2) = f.row(1);
	} else {
		flux_out.template topRows<2>() = f;
		flux_out.row(2) = f.row(0) * vPara;
	}
}

enum class WhichFlux {
	left,
	leftStar,
	rightStar,
	right
};

template<
	typename F_out,
	typename RawFlux, typename RawQ,
	typename Dl, typename Dr,
	typename Vnl, typename Vnr,
	typename Vpl, typename Vpr,
	typename Sl, typename Sr>
void setFluxComponents(
	F_out& f_out,
	const RawFlux& rawFluxL,
	const RawFlux& rawFluxR,
	const RawQ& rawQL,
	const RawQ& rawQR,
	const Dl& dL,
	const Dr& dR,
	const Vnl& vNormalL,
	const Vnr& vNormalR,
	const Sl& sL,
	const Sr& sR
){
	f_out = repl<2>(sL < 0 && sR > 0).select(
		repl<2>( useFluxLeft(dL, dR, vNormalL, vNormalR, sL, sR) ).select(
			hllFlux(
				sL, sR,
				rawFluxL, rawFluxR,
				rawQL, rawQR
			),
			hllFlux(
				sR, sL,
				rawFluxR, rawFluxL,
				rawQR, rawQL
			)
		),
		repl<2>(sL >= 0).select( rawFluxL, rawFluxR )
	);
}


template<
	typename Flux_out, Index bufCols,
	typename Dl, typename Dr,
	typename Vnl, typename Vnr,
	typename Vpl, typename Vpr,
	typename Sl, typename Sr,
	typename N>
void setFlux(
	Flux_out& flux_out,
	ArrayXNs<bufCols>& buf,
	const Dl& dL,
	const Dr& dR,
	const Vnl& vNormalL,
	const Vnr& vNormalR,
	const Vpl& vParaL,
	const Vpr& vParaR,
	const Sl& sL,
	const Sr& sR,
	const N& normal
){
	Index nCols { flux_out.cols() };
	assert( buf.cols() >= 2*5 );
	assert( buf.rows() >= nCols );

	Array2XsMap
		rawFluxL { buf.col(0).data(), 2, nCols },
		rawFluxR { buf.col(2).data(), 2, nCols },
		rawQL    { buf.col(4).data(), 2, nCols },
		rawQR    { buf.col(6).data(), 2, nCols },
		f        { buf.col(8).data(), 2, nCols };

	setRawFlux( rawFluxL, dL, vNormalL );
	setRawFlux( rawFluxR, dR, vNormalR );
	setRawQ( rawQL, rawFluxL, dL );
	setRawQ( rawQR, rawFluxR, dR );

	setFluxComponents(
		f, rawFluxL, rawFluxR, rawQL, rawQR, dL, dR, vNormalL, vNormalR, sL, sR
	);
	setFlux( flux_out, f, vPara, normal );

}


} // namespace hms::swe::hllc
